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Presenting simple coarse-grained models of isotropic solids and fluids in d = 1, 2 and 3 dimensions 
we investigate the correlations of the instantaneous pressure and its ideal and excess contributions 
at either imposed pressure (NPT-ensemble, A = 0) or volume (NVT-ensemble, A = 1) and for more 
general values of the dimensionless parameter A characterizing the const ant-volume constraint. The 
stress fluctuation representation of the compression modulus K in the NVT-ensemble is 

derived directly (without a microscopic displacement field) using the well-known thermodynamic 
transformation rules between conjugated ensembles. The transform is made manifest by computing 
the Rowlinson functional J^row also in the NPT-ensemble where = Kfo{x) with x = PidlK 

being a scaling variable, Pid the ideal pressure and foix) = x(2—x) a universal function. By gradually 
increasing A by means of an external spring potential, the crossover between both classical ensemble 
limits is monitored. This demonstrates, e.g., the lever rule Prow|x = A" [A -I- (1 — A)/o( 2 ;)]. 

PACS numbers: 05.20.Gg,05.70.-a,61.20. Ja,65.20.-w 


I. INTRODUCTION 


Strain and stress fluctuations. Among the fundamen¬ 
tal properties of any equilibrium system are its (general¬ 
ized) elastic constants characterizing the fluctuations of 
its extensive and/or conjugated intensive variables [H-Q- 
For instance, for an isotropic solid or fluid the volume 
and density fluctuations are set by the isothermal com¬ 
pression modulus K defined as 0 


K = V 



dP 

dP 

dV^ 

III 

T 


( 1 ) 


with V being the volume, N the particle number, p = 
N/V the particle density, F{T, V) the free energy, P the 
(mean) pressure and T the (mean) temperature. Albeit 
K may in principle be measured by fitting the pressure 
isotherm P{p, T) llOll. it is from the computational point 
of view important |lll. [l^ that this modulus may be ob¬ 
tained from the volume fluctuations at constant pres¬ 
sure (NPT-ensemble) and the pressure fluctuations at 
constant volume (NVT-ensemble) evaluated at the same 
state point, i.e. at the same mean temperature, density 
and pressure Q. In the NPT-ensemble K is obtained 
from the fluctuations SV = V — V of the instantaneous 
volume V around its mean value V = (V) using 


K = JVoilo = kBTV/ {5V^) 


( 2 ) 
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with fee being Boltzmann’s constant and where |o indi¬ 
cates that the average is obtained in the NPT-ensemble 
(A = 0 with A defined below). Equivalently, K may be 
obtained in a canonical NVT-ensemble using the stress 
fluctuation formula [mo, Od 


^ — ‘^Row'll — P ^Bornli ^F,ex|]^ (3) 

with |i indicating the NVT-ensemble (A = 1) and other 
definitions given immediately below. This formula was 
first stated for liquids in the 1950s by Rowlinson Q 
and later implicitly rediscovered by Squire, Hold and 
Hoover Q formulating the stress fluctuation formalism 
for anisotropic solids as summarized in appendix [Cl 
Affine contribution. The second term ryBorn in eq. @ 
represents the Born approximation [l^ for the interac¬ 
tion energy implied by an imposed infinitesimal strain as¬ 
suming affine microscopic particle displacements [1,0 0 ■ 
As reminded in appendix for pairwise additive poten¬ 
tials this “Born-Lame coefficient” becomes [HI 

?7Born = -P Hn)^ with h{r) = r (r ufr))' (4) 

and d being the spatial dimension, I an index labeling 
the interactions, r; the distance between two interacting 
particles, u{r) the pair potential and a prime denoting a 
derivative with respect to the indicated variable. As sug¬ 
gested by eq. Q, ’7Born is sometimes also called “hyper- 
virial” [n| . We note en passant that ryBom is a moment 
of the second derivative of u{r) and some care is required 
if JyBorn is Computed using a truncated potential [l5|, [l6j . 

Non-affne contribution. In general the Born approx¬ 
imation overpredicts the free-energy change. The over- 
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prediction is “corrected” by the stress fluctuation term 


VF,e 


= PV(SP', 


> 0 


( 5 ) 


with P = 1/k^T being the inverse temperature and Pex 
the instantaneous excess pressure which for pairwise ad¬ 
ditive potentials is given by Kirkwood’s virial m 

= with fi = -u'(ri) (6) 


being the central force between two interacting parti¬ 
cles. (Although Pex is used here in Rowlinson’s formula, 
eq. ©, at constant volume V = V, we have written it 
in a slightly more general form which is necessary if vol¬ 
ume fluctuations are allowed.) In particular from ref. @ 
it has become clear that stress fluctuation corrections, 
such as ?7F,ex|]^, do not necessarily vanish for T ^ 0. As 
we shall also illustrate in the present paper, this is due to 
the fact that the particle displacements need not follow 
an imposed macroscopic strain affinely SHMIl. How 
important the non-afhne motions are, depends on the 
system under consideration [2^. While the elastic prop¬ 
erties of crystals with one atom per unit cell are given 
by the Born term only, stress fluctuations are significant 
for crystals with more complex unit cells Q. They be¬ 
come pronounced for polymer-like soft materials and 
amorphous solids [1, d, ll7l - [^ . 

Fluctuations in conjugated ensembles. Focusing on 
the compression modulus K we emphasize in this report 
that the numerically more convenient stress fluctuation 
formalism may be obtained directly using the well-known 
thermodynamic transformation rules between conjugated 
ensembles [2^ . This point is crucial if the formalism 

is used in situations where no meaningful microscopic 
displacement field can be defined [id, [2^. Computing 
Rowlinson’s Prow for NPT-ensembles the general trans¬ 
form behind the formalism can be made manifest. Elab¬ 
orating a short comment [^ . we show that 

ProwIq = P ~\~ r^Bornlp ^F.exlp 

= K fo{x) with. fo{x) = x{2 - x) (7) 

being the universal scaling function, x = Pid/K the scal¬ 
ing variable and Pid the ideal pressure contribution. 

Generalized X-ensembles. It is straightforward to in¬ 
terpolate between the NPT- and the NVT-ensemble by 
imposing an external spring potential 

CextCfo) = ^(V- Pext) ' with Pext = C^ext (8) 

being the associated compression modulus introduced for 
convenience [s^. (Our approach is conceptually simi¬ 
lar to the so-called “Gaussian ensembl e” prop osed some 
years ago by Hetherington and others |3ll . l32j | generaliz¬ 
ing the Boltzmann weight of the canonical ensemble by 
an exponential factor Uext (E) oc of the instantaneous 
energy E.) Throughout this work it is assumed that 


AText > 0, i.e. Uext{V) reduces the volume fluctuations 
dll. Chosing the reference volume t4xt equal to the aver¬ 
age volume V of the isobaric system at imposed P allows, 
for symmetry reasons, to work at constant mean pressure 
irrespective of the strength of the external potential [sj . 
The volume fluctuations may then be characterized us¬ 
ing the dimensionless parameter A = Kext/{K + AText)- 
(Since AText > 0, we have 0 < A < 1 in the current study.) 
NPT-ensemble statistics is expected if the external po¬ 
tential does not constrain the volume fluctuations, i.e. 
A —>■ 0, while NVT-statistics should become relevant in 
the opposite limit for AText —>■ oo and A —>• 1. We shall 
monitor various properties, such as the Rowlinson for¬ 
mula pRowl;),, as a function of A and x = Pid/AT. We 
demonstrate, e.g., the simple lever rule 

-^RowU =7^[A+(l-A)/o(cr)] (9) 

which generalizes eq. m and eq. © to arbitrary A. 

Outline. In sect. m we consider theoretically vari¬ 
ous correlation functions of normal stress contributions 
in different A-ensembles. We begin by summarizing in 
sect. IB Al the transformation relations for fluctuations be¬ 
tween NVT- and NPT-ensembles. Equation © is derived 
in sect. ttmi We turn then in sect. IB El to the transfor¬ 
mation relations for general A and demonstrate eq. m 
in sect. IB El Our Monte Carlo (MC) and molecular dy¬ 
namics (MD) simulations of several simple coarse-grained 
models in d = 1, 2 and 3 dimensions are described in 
sect. imi Our theoretical predictions are then checked 
numerically in sect. lIVI Several well-known but scattered 
theoretical statements are gathered in the appendix. 


II. THEORETICAL CONSIDERATIONS 


A. Fluctuations in NVT- and NPT-ensembles 


As discussed in the literature [3, [Tl|, [2^ , a simple av¬ 
erage A = {A) of an observable A does not depend on 
the chosen ensemble, at least not if the system is large 
enough. (We do thus not indicate normally in which 
ensemble the average has been taken.) However, a cor¬ 
relation function (SASB) of two observables A and B 
depends on whether V or P are imposed. As shown first 
by Lebowitz et al. in 1967 [l^, one verifies that [13, 




dASBj 


K dAdB 


( 10 ) 


to leading order. We note that the left hand-side of 
eq. (Uni) must vanish if at least one of the observables 
is a function of B. In this case we have 


(^SASBj 


K dAdB 

Jv 


( 11 ) 


One verifies ior A = B = V that eq. ©1) is consistent 
with eq. @. For A = V and B = P one obtains imme- 
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diately the well-known relation Q 


-P (SVSP) 


= 1 . 


0 


Similarly, one obtains ioi A = B = V" that 

2 


(5(1/")' 


K d(V^} 


/3K 


0 ,3V y dP 

where the steepest-descent approximation 

\ 1/n 


j^2y2n-l 




V)=V 


( 12 ) 


(13) 


(14) 


for simple averages has been made and VjK = —dVjdP 
is used again. For the fluctuations of the inverse vol¬ 
ume l/y the latter result (n = —1) may be rewritten 
compactly using eq. as 



where we have changed « to the equal sign for large 
systems. That eq. (USD and eq. (USD become exact for 
F —> oo can be also seen by using that the distribution 
of V in the NPT-ensemble is Gaussian. With A = F" 
and B = P one gets similarly 


<Jd(F")(5Pj>|^ 


K a(F") 
W dP 


-nV^-^ksT 


(16) 


to leading order for F —>■ oo using the same approxima¬ 
tion as above. With eq. m this gives for n = — 1 the 
convenient cumulant 


= 1 . 


ysvsp 


(17) 


The cumulants eqs. (I15I17I) have been used in the com¬ 
putational part of our work to check the precision of the 
barostat and to verify whether our configurations are suf¬ 
ficiently large for the investigated state point. 


B. Transformation of pressure auto-correlations 


Returning to eq. (USD, this implies tor A = B = P the 
transformation of the pressure fluctuations 


f3v ysp'^ 


^ = 13V (SP'^ 


-K, 


(18) 


i.e. K may be obtained by measuring the pressure fluc¬ 
tuations in both ensembles. As we shall show in para¬ 
graph III Cl the numerically more convenient Rowlinson 
expression ProwIi for K can be derived directly from 


eq. (fT51) [l3l| . In the following we use the more concise no¬ 
tation rjp = (3V{6P^) for the pressure fluctuations. 77 f|o 
is also called the “affine dilatational elasticity” tja [I3| for 
reasons which will become obvious in sect. Ill Cl (See also 
appendix [All Since AT > 0 for a stable system [ 3 , eq. (fT5l) 
implies rjA = ?7fIo > ^fIi- Depending on the disorder, 
77 f|]^ is, however, not a negligible contribution. To see 
this let us remind that K can also be determined from 
the (F, P)-data measured in an NPT-ensemble using the 
linear regression relation 


if= Preglo= -V{5V5P) 



(19) 


Please note that using eq. (IT^ this reduces to eq. ©■ 
Associated to Preglg is the dimensionless regression coef¬ 
ficient 


'"I'eglo ~ 


{svsp'y 


0 


/ 



( 20 ) 


which can be also further simplified using eq. ((HD- In¬ 
terestingly, using eq. @ and eq. (ITSl) one sees that 


Creglo = VK/va = \J3- »7fIi Ma, (21) 

i.e. the regression coefficient obtained at constant pres¬ 
sure determines the pressure fluctuations at constant vol¬ 
ume. Only if the measured (F, P) are perfectly corre¬ 
lated, i.e. Creglo = this implies K = rjA and r]p\.^ = 0. 
In fact, for all non-trivial systems one always has 


Creglo K < rjA and 77 fIi > 0, (22) 

i.e. the affine dilatational elasticity rjA sets only an upper 
bound to the compression modulus and a theory which 
only contains the affine response must overpredict K. 


C. Rowlinson’s formula rederived 


MC-gauge. There is a considerable freedom for defin¬ 
ing the instantaneous value of the pressure P = Pd -f Pex 
as long as its average P = Pid -f Pex does not change [111. 
It is convenient for the subsequent derivations and for 
our MC simulations (and not in conflict with the also 
presented MD simulations) to define the instantaneous 
ideal pressure Pd by 


Pid=kBTN/V (MC-gauge). 


(23) 


Within this “MC-gauge” the thermal momentum fluctu¬ 
ations are assumed to be integrated out. This leads to 
the usual prefactor 


F'^ = exp 


-f3{-kBTNlog{V)) 


(24) 


of the remaining partition function. The effective Hamil¬ 
tonian 'Hsiy) of a state s of the system thus reads 

UsiV) = -kBTNlogiV) + Us{V) + const (25) 
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where the first term on the right hand-side refers to the 
integrated out momenta and the second to the total ex¬ 
cess potential energy Us{V) expressed as a function of 
the instantaneous volume as shown in appendix lAl 
Non-ajfine contribution. An immediate consequence 
of the MC-gauge is, of course, that the fluctuations of 
Pid vanish for the NVT-ensemble and that, hence, 


iiY\, = mM^ = pv 


(26) 


with the instantaneous excess pressure Pex being com¬ 
puted using Kirkwood’s expression, eq. (I6|). According 
to eq. m the correlation coefficient Creg is thus a func¬ 
tion of ?7F,exli /VA and vice versa. 

Affine (Born) contribution. The task is now to com¬ 
pute the pressure fluctuation in the NPT-ensemble. We 
note first for the NPT-weight of a configuration at vol¬ 
ume V that 


'^-0{-Hs{V)+PV) 


ffiP _ P)e-P(f^AV)+PV) (27) 


where P = (P-gauge) (28) 


defines the instantaneous total pressure P = Pid -I- Pex- 
Note that this definition is consistent with the MC-gauge, 
eq. (l23)) . for Pid and eq. (ESI). As shown in appendix |B1 
Pex = —U({V) is also consistent with the Kirkwood ex¬ 
cess pressure, eq. (0), for pair potentials. Using eq. EZD 
the second moment ((P — P)^) can be readily obtained 
by integration by parts. This yields 




u(p"(U))|o 


(29) 


where apriori the average is understood to be taken over 
all states s of the system and all volumes V at imposed P. 
(The boundary terms for the integration by parts over V 
can be neglected for sufficiently large systems since the 
NPT-weight for V gets strongly peaked around V.) It 
is of importance that the fluctuation ?7 f|o has thus been 
reduced to a simple average. This allows its computation 
more conveniently by NVT-ensemble simulations. Using 
eq. EH) it follows further that rjx = pa, id + ??A,ex with 


^A.id 


^A,ex 


u(-fcBTVlog(U)")[ 


Pid 


(30) 

(31) 


where we have used for the ideal contribution that for suf¬ 
ficiently large systems U(l/U^) « (1/U) ~ 1/U. As seen 
from the affine excess energy discussed in appendix it 
follows for pair potential interactions that 


^A,ex — ^Born 4“ Pex; (^^) 

i.e. both coefficients ryA.ex and ?7Born are equivalent. We 
stress again that Pid, Pex, VBom and r]A,ex are simple 
averages and can thus be evaluated readily in both en¬ 
sembles using eq. dn. Substituting Gcjs. into 

the transform eq. this finally confirms eq. ED- 


D. Correlations at constant P 


We focus now on the fluctuations of the pressure con¬ 
tributions in the NPT-ensemble. According to eq. (l32ll 
we have P -I- ?7Born = ^a = ?7fIo- If the Rowlinson func¬ 
tional pRow is measured at imposed P this implies 


ProwIq — ^fIo VP,e 

= pv (dAd 


2/3U (<5Pid5Pe: 


.(33) 


As a next step we demonstrate the relations 


^7F,idlo=/3u(^Pid)[ 

»?F,mix|o = Py (^SPidSPe^^ ^ 


K (34) 

Kxil-x) (35) 


(with X = Pid/K being again the scaling variable) 
from which eq. ED is then obtained by substitution into 
eq. (l33)l . Remembering eq. (ESD, eq. (IMl) is obtained from 
the fluctuations of the inverse volume, eq. (IT^ . The rela¬ 
tion eq. (EID describing the coupling of ideal and excess 
pressure is implied by eq. EID for n = — 1 and using 
P = Pid + Pex and eq. dMD ■ We emphasize that eq. ED or 
eq. (13311 do not completely vanish for finite T, i.e. finite 
X, as does the corresponding stress fluctuation expression 
for the shear modulus G at imposed shear stress r (T^ . 

Having thus demonstrated eq. (IMl) and eq. (1351) and 
using the already stated relation for the total pressure, 
eq. (HHD, one confirms finally for the fluctuations of the 
excess pressure Pex that 


flF.exlo = Py {SPL) 


= T?A.ex - K x{l 

0 


x) 


(36) 


with x = Pid/AT being the reduced ideal pressure. 


E. Fluctuations in different A-ensembles 

Introduction. The transformation relations for sim¬ 
ple averages and fluctuations between the standard con¬ 
jugated ensembles have been given in the 1960s by 
Lebowitz, Percus and Verlet [2^. Focusing on the volume 
V as the only relevant extensive variable and the conju¬ 
gated (reduced) pressure fP of the system as the only 
intensive variable [13, we rederive now their results 
for generalized A-ensembles. We illustrate several points 
made by the histograms presented in fig. [T] for simple ID 
nets of harmonic springs as described in sect. IIIII below. 

Histogram. Let us assume a normalized distribution 
p{V) with a sharp and symmetric maximum at the mean 
volume V. Examples for such histograms are given in 
panel (a) of fig. [T]for several values A = + K). 

The standard NPT-ensemble corresponds to the value 
A = 0. Note that the volume V of the NPT-ensemble at 
P = 0 is taken as the reference volume 14xt of the exter¬ 
nal spring potential. Due to this choice neither the mean 
volume V nor the pressure P do change with increasing 
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FIG. 1: Simple ID nets without noise {5k = 0) at tempera¬ 
ture T — 0.01 and pressure P — 0 (mean volume V = 10^, 
number density p = 1) for different values A of the exter¬ 
nal spring as indicated: (a) Normalized histograms p{x) of 
the rescaled instantaneous volume x = V/V — 1. The dis¬ 
tributions are Gaussian as shown by the line indicated for 
the standard NPT-ensemble (A = 0). (b) Pre-averaged pres¬ 
sure SP{V) traced as a function of the rescaled instantaneous 
volume X. The bold line shows eq. dSZl) for the compression 
modulus K = 1 of the system. Inset: The mean-squared vol¬ 
ume fluctuations SV^ decrease linearly with A according to 
eq. (gH) or eq. (1481) as shown by the solid line. 


A, i.e. all ensembles correspond to the same thermody¬ 
namic state. All distributions are Gaussian, becoming 
sharper with increasing A. (The limit A —>■ 1 corresponds 
to a Dirac 5-function.) The width of the distribution 
around its maximum being characterized by the param¬ 
eter A, we introduce the notation 5V^ = {SV"^) for the 
mean-squared width. The width 5V^ of simple ID nets 
is presented in the inset of panel (b) of fig. [T] 

Observables. We write A{V) for the expectation value 
of an observable A at a given state s of the (full) phase 
space of a system at volume V. Let us first focus on 
the average A{V) of A{V) taken over all states s. As an 
example we show in panel (b) of fig. [T] the mean pressure 
P{V) averaged over all systems s of a A-ensemble found 
at a given instantaneous volume V. As one expects, P{V) 


is seen to decrease linearly to leading order, 

5P{V) = P{V) - P{V) « SV = V-V, (37) 

as stressed by the bold line. Importantly, P(V) does not 
depend on A. (As one expects, the statistics deteriorates 
for {V — V)'^j6V^ ^ 1.) More generally, a Taylor expan¬ 
sion of A{V) around V yields 

A{V) = Up -|- uiSV —a2bV'^ -I-... (38) 

where a coefficient a„ denotes the nth derivative of A{V) 
with respect to V taken at V = V. 

Simple averages. The average taken over all properly 
weighted volumes of the investigated A-ensemble thus be¬ 
comes to leading order 

(a( 1 /)) = J dVp{V)A{V) « ap + (39) 

with Op = A{V). We have used here that p{V) is sym¬ 
metric around V which implies that the linear term in 
eq. (|55)) must drop out. The difference of a simple av¬ 
erage taken at A < 1 and taken in the NVT-ensemble 
(A = I) is then given by the second derivative 02 times 
the mean-squared width 5V^ which increases at most lin¬ 
early with y for A = 0. Since 02 does not depend on 
the ensemble, the correction must decay at least as fas t 
as for the NPT-ensemble discussed in the literature [^ . 
Hence, simple averages computed for any ensemble with 
0 < A < 1 become rapidly indistinguishable. 

Averaged fluetuations. As a next step let us consider 
the correlation function {SA6B) = (AB) — {A){B) of the 
observables A{V) and B{V). Using again the symmetry 
of the distribution p{V) around V it is seen that 

(SASB) = (gpbp) - (ae)( 6 p) -I- aibiSV^ (40) 

where the coefficients a„ and stand for the respective 
derivatives of A{V) and B{V) taken at V = V used for 
the Taylor expansion around V. The underlined term 
vanishes since ap = A{V) and 6 p = B{V) are constants. 
Note also that {6A6B) —>■ 0 for A —>■ I. This is due to 
the fact that we have correlated here the pre-averaged 
observables A{V) and B{V) instead of the expectation 
values A and B which depend not only on the volume 
V but also on the state s of the system. Replacing thus 
A ^ A, B ^ B, On ^ dn and —)> bn we have thus 
in addition to average over all possible states and the 
underlined term in eq. (HOI) remains thus finite 

(SASB) = {6do5bo) + {dibi)6V^. (41) 


Since (ai 6 i) is a simple average, this result can be refor¬ 
mulated using a more natural notation as 


SASBj 


= (^SASB'j 


dAdB 
dVdV ^ 


(42) 
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where |a indicates that the average is taken over properly 
weighted volumes in a general A-ensemble and 11 indicates 
the NVT-average (A = 1) taken at the maximum of the 
distribution p(y). Using eq. for the NPT-limit (A = 
0 ) this is seen to be identical to eq. (2.11) given in ref. (^ . 
Substracting this reference from the general A case yields 


(^SASb'^ ^ = (^SASB'j 


/ ^ 9 ^ 9x dA dB , , 

„ - m - sva (43) 


Equation (1431) can be further simplified using 

V V VX 

Jk ~ P{K + i^ext) “ W 


dU? - 5V^ = 


(44) 


for the difference of the volume fluctuations in both en¬ 
sembles. We have used here eq. (HJ for the NPT-ensemble 
and the fact that for general A the external spring is par¬ 
allel to the system, i.e. the effective modulus must be the 
sum of the system modulus K and spring modulus if ext ■ 
That eq. (l4^ holds is confirmed by the data presented in 
the inset of panel (b) of fig. [TJ Following ref. we also 
rewrite the derivatives with respect to the volume V as 
derivatives with respect to the pressure P of the system 


The latter result is consistent with eq. m which is thus 
shown to hold for all A < 1. For A = V and .B = P it is 
seen that 


-p (^svsP'j 


-p (^svsP) ^ 


KXdV_ 
~dP ~ 


(51) 


generalizing thus eq. m- More generally, one sees for 
A = U" and B = P that 


-p (^diV^)SP^ 


= (1 - A) 


(52) 


as expected from eq. (HU). One confirms using eq. dUD 
that eq. (II3 must hold for all A < 1. Interestingly, the 
already mentioned linear regression formula Preg is seen 
using eq. (H51) and eq. m to become 


-V (SVSP) 


/(sr^ 


k^TV 

5Vo^ 


= K 


. 

independent of the ensemble used. As an alternative to 
the strain fluctuation relation eq. ((l^ . this allows thus 
the determination of K for all A < 1. The dimension¬ 
less correlation coefficient Creg associated to p-eg depends 
however on A 


dAdB _ /'dPV dAdB _ /kV dAdB 
dVdV~[dv) 'dPdP~[y) Pp'dP 


where we have used dPjdV = —KjV in the last step. 
Using the above three equations this yields finally 


{5A5B) 


{5A5B) 


KX dAdB 
W dPdP' 


(46) 


which compares the correlations in a general A-ensemble 
(0 < A < 1) with the correlations in an NPT-ensemble 
(A = 0). Note that for A = 1 this is consistent with the 
original transformation, eq. (HOD, derived in ref. [^ . 


F. Correlations for generalized A-ensembles 

Using eq. (I46D we restate first several correlations given 
above where at least one of the observables A and P 
is a function of the instantaneous volume V. Since for 
A < 1 volume fluctuations are not (completely) sup¬ 
pressed, eq. (HU cannot be generalized. Instead the al¬ 
ready demonstrated results for A = 0 or A = I are used. 
For A = P — V this yields, e.g., 

{svA 


restating thus eq. (H^ . More generally, one confirms for 
i = P = U" that 


{5V^) 


KX dV 


0 PV \dP 


(I - A) 
kBTV 
K + Pext 


(47) 

(48) 

(49) 


I - /Z / 

V V 1 - 


(54) 


where we have used eq. dSSD demonstrated below. Note 
that P/?7A = 1 implies (as before) Cregl;^ = 1 for all 
A < 1. For Ar/?7A < 1 the correlation coefficient decreases 
continuously from its maximum \JK/rj^ at A = 0 to 
zero for A = I. As one would expect, this shows that the 
more the volume fluctuations are suppressed by the exter¬ 
nal constraint, the more V and P must decorrelate. For 
the transformation of the total pressure auto-correlation 
r]F = PV{6P^), eq. (I46D simply implies that 


VfIx = ^fIo - pa = 77a - 


PPext 

K + Pext ■ 


(55) 


Since rjA = Pid + ?7A,ex is a simple average computable 
in any ensemble, eq. (|55D may be also used for the de¬ 
termination of K. Using again the MC-gauge and that 
P = Pd + Pex, it is seen from eq. (1551) that the Rowlinson 
functional ProwIa = P + ?7 Born — ?7F,ex|;^ should become 


pRowl;^ — KX 77F,id|;^ 4” 2 77F,mix|;^ (^^) 

generalizing thus the Rowlinson formula, eq. (HD , for A = 
1 (remembering that 77F,id and ryF.mix must vanish in the 
NVT-ensemble) and eq. (1551) for A = 0. Considering as 
a next step the fluctuations of the ideal pressure 77F,id in 
the MC-gauge, eq. (1501) implies that 


77F.idU = Py (<^Pd) = P P (1 - A). (57) 


Since P = Pid+Pex the latter result together with eq. (1551) 
allows to generalize eq. (1581) for the correlations between 
ideal and excess pressure contributions. This shows that 


(d(U")2) 


= n^V 


2-\/-2n — 2 


(1-A). 


(50) 


= Py {^P<iPe.)\^ = K x{l-x) (1-A). (58) 
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Our central result eq. m announced in the Introduc¬ 
tion is thus simply obtained by substituting eq. and 
eq. (1551) into eq. (1551) . Finally, using eq. (|55|) together 
with eq. (|57|) and eq. (1551) one verifies for the fluctuations 
of the excess pressure fluctuations that 

^F,ex|^ = VA -^RowI;, (^1^) 

as obvious from eq. (0. The latter result reduces to 
the Rowlinson expression eq. for A = 1 and using 
Va = ?7A,ex + Kx to eq. (|55)) for A = 0. 


III. SOME ALGORITHMIC DETAILS 

Introduction. In order to check our predictions we 
sampled by MC and MD simulation [IJ, [T^l various 
model systems for solids and glass-forming liquids in 
d = 1, 2 and 3 dimensions. Periodic boundary conditions 
are used and all systems are first kep t at constant pres¬ 
sure P using standard barostats [ll| as specified below. 
After equilibrating and sampling in the NPT-ensemble, 
the volume fluctuations are suppressed either by impos¬ 
ing V = V or hy means of a finite spring potential. Var¬ 
ious simple averages, such as the Born-Lame coefficient 
??Born, and fluctuations, such as the excess pressure fluc¬ 
tuation r?F,ex or the Rowlinson functional Jfiiow, are com¬ 
pared for all available A. As expected, all simple averages 
are found within numerical accuracy to be identical. As 
discussed in sect. lYl fluctuations are found to transform 
following the predictions based on eq. (nni) and eq. (|551). 
Specifically, we have verified that eqs. (I15I17I5II) hold 
to high precision for all A. About N = 10^ particles 
are typically used. Lennard-Jones (LJ) units are used 
throughout this work and is set to unity. 

One-dimensional spring model. As already seen in 
fig. [TJ the bulk of the presented numerical results has 
been obtained by MC simulation of permanent nets of 
ideal harmonic springs in strictly d = 1 dimension [3^ . 
We use a potential energy 

^ = ( 60 ) 
I 

with xi being the distance between the connected parti¬ 
cles. The reference length Ri of the springs is assumed 
to be constant, Ri = R = 1, and the spring constants ki 
are taken randomly from a uniform distribution of half¬ 
width Sk centered around a mean value also set to unity. 
We note for later reference that this implies 


{1/ki) \og{{l + 6k)/{l-6k)) ^ ' 

as indicated by the bold line in fig. [2j Only simple net¬ 
works are presented here where two particles i — 1 and i 
along the chain are connected by one spring I = i, i.e. all 
forces fi along the chain are on average identical. 


Glass-forming liquids. Our two-dimensional (2D) sys¬ 
tems are polydisperse Lennard-Jones (pLJ) beads as de¬ 
scribed in ref. [I3l | . The reported three-dimensional (3D) 
systems refer to binary Kob-Andersen (KA) mixtures 
sampled by means of MD simulations taking advantage 
of the LAMMPS implementation [s^. Starting from the 
liquid limit well above the glass transition temperature 
Tg, both system classes have been quenched [l^ deep 
into the glassy state at very low temperatures T <C Tg as 
may be seen from fig. [3] [1^ . 

Barostats. The results reported for d = 1 and d = 2 
have been obtained using local MC moves for the parti¬ 
cles and global MC moves for the barostat [T^. As de¬ 
scribed elsewhere [ll|) [a, an attempted volume change 
5V = Vnew — I4id is accepted if ^ < exp(—/3<5G) with ^ 
denoting a uniformly distributed random variable with 
0 < ^ < 1 and 

5G = 5E + SVP - kBTN log(I4ew/I4id) ■ (62) 

The first term 6E stands for the energy difference asso¬ 
ciated with the affine displacements of all particles, the 
second term 6VP imposes the pressure P and the log¬ 
arithmic contribution corresponds to the change of the 
translational entropy, i.e. the change of the integrated 
out momentum contribution discussed above, eq. (1241) . 
While a broad range of pressures has been sampled for 
the ID nets, the 2D pLJ beads have been kept at only one 
pressure, P = 2, for which a glass-transition temperature 
Tg « 0.26 has been determined 

For more general A-ensembles the increment SUext(,V) 
of the external spring defined in eq. ® is simply added to 
SG. We assume throughout this work that I4xt = V for 
the reference volume of the external spring with V being 
the mean volume of an NPT-ensemble of a given pressure 
P. Due to the symmetry of the fluctuations around V 
for the system and the external spring, this is sufficient 
to keep this pressure constant for all A [s^. 

For our MD simulations of the KA model we have used 
the Nose-Hoover barostat (“fix npt command”) provided 
by the LAMMPS code [s^. Following Kob and Andersen 
[38l| a constant pressure P = I has been imposed for all 
temperatures. This choice corresponds to glass-transition 
temperature Tg « 0.4 [T^ 


IV. COMPUTATIONAL RESULTS 

As already noted above, various properties have been 
compared for different boundary conditions as character¬ 
ized by the parameter A while keeping the system at the 
same state point. We focus first on the comparison of 
NPT- (A = 0) and NVT-ensembles (A = 1) before we 
turn in sect. lIVEl to the more general A-ensembles. 









FIG. 2: Compression modulus K for ID nets at temperature 
T = 0.01 and pressure P — 0 computed using the rescaled vol¬ 
ume fluctuations TvoiIq (filled spheres), the Rowlinson stress 
fluctuation formula (crosses), the difference between 

the total pressure fluctuations in both ensembles (squares) 
and the fluctuations of the inverse volume F)d/ ’iF.idIo (large 
spheres). The compression modulus decreases strongly with 
(5fc (solid line). Also indicated are the “afhne” contribution 
^f|q = l-|-2Pid to the compression modulus, eq. (IMJ, and the 
“non-afhne” correction which is seen to increase with Sk 
according to eq. (1671) as indicated by the dash-dotted line. 



FIG. 3: Compression modulus K vs. temperature T compar¬ 
ing again the values obtained using PvoiIq, Prow|j, eq. (fT^ 
and eq. (1341) . The upper data refer to the results obtained 
for the KA model in 3D [l^ . [3^ which have been shifted up¬ 
wards (factor 4) for clarity, the data in the middle to systems 
of glass-forming 2D pLJ beads at P = 2 [a and the lower 
data to simple ID nets of harmonic springs at P = 0 with 
identical spring constants (Sk = 0). 


A. Compression modulus K 

Comparing different methods. As shown in figs. ^ 
m the compression modulus K may be determined ei¬ 
ther using the volume fluctuations in the NPT-ensemble, 
eq. dSl), or using Rowlinson’s stress fluctuation formula, 
eq. ([ni), for the NVT-ensemble (crosses). The same val¬ 
ues of K are obtained from the transform of the pressure 



FIG. 4: Compression modulus K vs. pressure P for ID nets 
without noise {Sk = 0) at T = 0.01. The compression modu¬ 
lus decreases with P simply since the volume V of the system 
linearly decreases as F = Fref(l — P/Kmi) with Aref being 
the compression modulus at P = 0 where V = Fef. If K is 
rescaled with Vief/F (spheres) this yields a pressure indepen¬ 
dent material constant Aref. 

fluctuations, eq. (HU, and from the ideal pressure fluctu¬ 
ations ?7F,idlo: H34L as indicated by the large spheres 

which thus confirms both relations. 

Harmonic spring networks. As seen in fig.[21 the com¬ 
pression modulus of the ID nets decreases strongly with 
6k. To understand the scaling indicated by the bold line 
let us consider a chain of harmonic strings. Since the av¬ 
erage force acting on each spring is given by the imposed 
pressure P = —/ext/A, this implies an average length 
{xi) = R — PA/ki for a spring constant ki and an av¬ 
erage volume V — A J2i{xi). A pressure increment SP 
thus leads to a volume change 6V = —ASPJ2i ^/^i- 
the compression modulus K = —V6P/6V at pressure P 
this yields 

A = Aret - P with Aret = Pref/ (1/^i) (63) 

being the compression modulus of the unstressed refer¬ 
ence system at P = 0 and p^ef = l/AR the correspond¬ 
ing density. The compression modulus A = Aj-et at zero 
pressure is thus inversely proportional to the average in¬ 
verse spring constant. For our uniformally distributed 
spring constants eq. (1611) thus implies that A must van¬ 
ish as a cusp-singularity for dk ^ 1 [i^. Also indicated 
in fig. are the “affine” contribution pf|q to A (further 
discussed in sect. IIV B|) and the “non-affine” contribution 
PfIi which is seen to increase with Sk. The decrease of A 
is thus due to the increase of the non-afhne contribution. 

Temperature dependence. Figure [3] shows the temper¬ 
ature dependence of the compression modulus for our 
three model systems. As one would expect, A decreases 
with T for the bead systems kept at same pressure reach¬ 
ing the ideal gas compressibility A = Aid = Pid = P for 
large temperatures. As predicted by eq. (1531) for strictly 
harmonic spring chains, the compression modulus A is 
found T-independent for all ID networks and this irre- 
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spective on the values 5k and P sampled (not shown). 

Pressure dependenee. The pressure dependence of the 
compression modulus for ID nets is investigated in fig. 01 
The dashed line indicating the linear behavior predicted 
by eq. ([Ml) for chains of harmonic springs perfectly fits 
the measured data points. A comment is in order here: 
Since the springs are permanently fixed and ideal, the 
intrinsie mechanical properties of the systems do not 
change with the external load. In fact, the second deriva¬ 
tive of the free energy F{T, V) with respect to the volume 
V becomes constant for these highly idealized systems 


K _ d^F{T, V) 
V ^ 9D2 


= const 


ATref 

Vre{ 


(64) 


with V(.ef being the reference volume of the system at 
zero pressure. The difference between the thermody¬ 
namic compression modulus K and the constant ATj-ef 
simply arises since in all thermodynamic relations, such 
as eq. (Cl), the volume V of the current state is taken 
to make the modulus system-size independent (inten¬ 
sive) and not a reference volume at a certain pressure 
P. Since for an idealized elastic body as our ID nets 
V = V)ef(l — P/Kref), this implies the linear relation 
K = Kj-ei — P (dashed line). Note that for solids the ap¬ 
plied pressures are for once small compared to the moduli 
and the systems can be taken linearly and without (or at 
least with negligible) plastic rearrangement from a zero- 
pressure reference to a finite P. It is thus possible to 
define elastic material properties such as ATj-ef which are 
independent of the applied external stresses. If possible, 
this is a conceptionally and practically very useful sepa¬ 
ration of conditions and effects. However, it is obviously 
pointless to take the zero-pressure volume (or the vol¬ 
ume of any other pressure) as a reference for a colloidal 
or polymeric system from which other pressure states are 
described by a linear extrapolation. Please note that be¬ 
ing in principle not very deep (being a matter of defi¬ 
nitions), this issue creates significant confusion between 
people from different communities. See appendix [C] for 
addition comments on this issue. 


B. Affine compressibility contribution 

The inset of fig. [5] presents for our three model sys¬ 
tems the correlation function rjp = PV{6P^) of the total 
pressure in the NPT-ensemble measured as a function 
of temperature while keeping the mean pressure P con¬ 
stant. As shown in sect. m the “affine dilatational elas¬ 
ticity” TjA = ?7fIo yields the leading contribution to the 
compression modulus K. We note first that rjA increases 
with decreasing temperature for our two glass-forming 
liquids in d = 2 and d = 3 dimensions (due to the in¬ 
creasing repulsion of the LJ beads) levelling off below 
the glass-transition temperature Tg of each model. At 
contrast, tja is seen to increase monontonously with tem¬ 
perature. This qualitatively different behavior needs to 
be explained. We remind first that rjA can be expressed 



FIG. 5: Affine dilatational elasticity t^a determined from the 
pressure fluctuations ryrlQ in the NPT-ensemble. Inset: rjA as 
a function of T for ID nets {Sk = 0, P = 0), 2D pLJ beads 
(P — 2) and 3D KA binary mixtures (P = 1). Main panel: 
Data collapse for t^a -f P as a function of Pid(T) for ID nets 
for a broad range of parameters as indicated. The bold fine 
corresponds to the expected behavior for harmonic springs. 


by the simple average, eq. ([29]), which can be evaluated 
in any ensemble. For the pair potentials we focus on here 
this yields rjA = PeK+mom = Pid+r]A,ex with ryBom being 
the Born-Lame coefficient, eq. o, as we have explicitly 
checked for all models M- Using eq. (IA13I) one sees that 
for harmonic springs in d = 1 dimensions 

VA = Pid + h (xf) (65) 

for a given quenched realization of spring constants ki. 
Since in the NPT-ensemble each spring is uncorrelated, 
we can use that (xf) — {xiY = k-^T/ki for the thermal 
fluctuation of every spring 1. This implies that rjA = 
‘^Pid+J2i /^- Since the average length of a spring 

at constant pressure P is given by (a;/) = i? — PA/ki and 
since the mean volume V is the sum of averaged lengths 
(xi) times the surface A, it follows using eq. (l63l) for the 
compression modulus K that 


_OD , I,- r, , {ki){l/ki)-l 

77A - 2Pid + dfref P + ^ (66) 

with ATref = Prefd?^/(l/fci) and pref = 1 /AR characteriz¬ 
ing the zero-pressure reference. For systems of identical 
springs this reduces to ija = 2Pid + ATref — P- With 
ATi-ef = I this is shown by the bold line in the main panel 
of fig. [S] For polydisperse systems at zero pressure we 
have instead tja = 2Pid + PrefR^iki.) as already shown by 
the dashed line in fig. [2] Within the units used this cor¬ 
responds numerically also to pA = 1 + 2Pid. Hence, all 
ID nets without noise or without applied pressure should 
collapse on the same master curve if pA + A is plotted as 
a function of the ideal pressure Pid(T). This is confirmed 
by the data presented in the main panel of fig. |5l 
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FIG. 6: Characterization of non-afhne elastic response. Main 
panel: Excess pressure fluctuation as a function of 

Pid{T) for ID nets for different Sk and P. Systems with no 
{Sk = 0) or little polydispersity and all high-temperature sys¬ 
tems collapse on the dashed line indicating ryp.exlj = 2Pid. 
Inset: Dimensionless regression coefficient Creglo for all three 
model systems as a function of t^pIj^ / ’IpIq = ’Ip.exl^ /va con¬ 
firming eq. as indicated by the bold line. 


C. Non-afflne compressibility contribution 

The non-affine deviations from the affine Born con¬ 
tribution rjA to the compression modulus K are given 
according to the Rowlinson formula, eq. by the 
fluctuation ?7F,exli of the excess pressure in the NVT- 
ensemble. Focusing on one low temperature we have al¬ 
ready seen ?7F.exli in fig. [Has a function of 6k. Plotted as 
a function of the ideal pressure Pid(T'), the main panel 
of fig. [6] presents ?yF,ex|i for various 6k and P as indi¬ 
cated. Since ?yF,ex|i = VA — K, it follows from eq. (EH) 
and eq. (IS51) characterizing, respectively, the compression 
modulus and the affine dilatational elasticity, that for ID 
nets of harmonic springs 


^F.exl]^ 


2Pid + 


l-P{l/ki) 


(67) 


where we have used that the mean spring constant (fc/), 
the reference length R of the springs and the surface A 
are all arbitrarily set to unity. Assuming P = 0, eq. (EH) 
is traced in fig. [5] as a function of 6k (dash-dotted line). 
Since the second term vanishes for identical springs, this 
relation implies »7F,ex|]^ = 2Pid for all temperatures and 
pressures as shown in the main panel of fig. |6| by the 
dash-dotted power-law slope. This is confirmed by the 
data presented for 6k = 0. As one expects, this also 
gives the high-temperature limit for all systems. Note 
that r7F,ex|]^ only vanishes in the low-P limit of strictly 
identical springs. Since the second term in eq. EH) is a 
finite constant for 6k > 0, the corresponding data must 
thus level off for T —> 0. Using a simple example we 
have thus confirmed the more general finding by Lutsko 



FIG. 7: Stress fluctuations in the NPT-ensemble. Large 
spheres refer to 2D pLJ beads for P = 2, squares to 3D KA 
beads for P — 1, all other symbols to ID nets for different 
Sk and P as indicated. Main panel: y = ?7F,mix|o/Pid as a 
function of the reduced ideal pressure x = Pa/K confirming 
eq. (|35|. Inset: Similar scaling for the reduced correlation 
function y = (?7A,ex — ??F,ex|o)/Pid using the same symbols as 
in the main panel. The bold line represents eq. pel). 


Q that the non-affine contributions to the elastic moduli 
need not necessarily vanish in the low-T limit. Inter¬ 
estingly, due to the P-dependent denominator the low-T 
non-affine contribution becomes more and more relevant 
for larger pressures. It becomes thus increasingly diffi¬ 
cult to reach the asymptotic high-T limit (not shown) 
(4ll |. As stressed in sect. imi the correlation coefficient 
Creg of the instantaneous volumes and pressures in the 
NPT-ensemble is set by the reduced non-affine contribu¬ 
tion ?7 f|i / ?7f|o = ?7F,ex|]^ /rjA. This relation is verified in 
the inset of fig. [5] presenting for all three studied model 
system a perfect collapse of the data on the master curve 
(bold solid line) predicted by eq. (EH). Hence, the coef¬ 
ficient Creg may be seen as a dimensionless and properly 
normalized “order parameter” characterizing the relative 
importance of the non-affine displacements. This also 
shows that it is equivalent to specify for a system either 
(A:,77a), (r?A, m ,ex|x ) or (AT, Creg). 


D. Correlations in the NPT-ensemble 

As shown in fig. [71 we have checked the predicted 
correlations between the ideal and the excess pressure 
fluctuations jyF.mixlo the auto-correlations of the ex¬ 
cess pressure ?iF,ex|Q- To make all investigated mod¬ 
els comparable the reduced correlation functions y = 
ffF.mixlo /Pid (main panel) and y = (?7A.ex - ??F,exlo)/Pid 
(inset) are traced as a function of the reduced ideal pres¬ 
sure X = Pid/AT with K as determined independently 
above. Please note that compared to our theoretical pre¬ 
dictions we have scaled the vertical axes of the data not 
with AT, but using the ideal pressure Pm, i.e. the predic¬ 
tions have been divided by x. This was done for presen- 
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tational reasons since thus the low-T asymptote becomes 
a constant. A perfect data collapse on the (rescaled) pre¬ 
diction y = 1 — X (bold line) is observed for all systems. 
As shown in ref. [23 a similar data collapse is achieved 
for the Rowlinson formula if y = -^Row Iq /^id is plotted as 
a function of x. (Since a related scaling plot is presented 
in sect. II V El we do not repeat this figure here.) Interest¬ 
ingly, the latter scaling does not depend on the MC-gauge 
which has been used above to simplify the derivation of 
eq. 0. Please note that for our glass-forming liquids 
it is not possible (for the imposed pressures) to increase 
X beyond unity {K > Pid) and the deviations from the 
predicted plateau y = fo{x)/x ~ 2 for a: 1 are thus 

necessarily small and can be easily overlooked (T^ . 


E. Generalized A-ensembles 

Having characterized the behavior in the NPT- 
ensemble (A = 0) and the NVT-ensemble (A = 1) we 
turn now to the general A-ensembles. As described at the 
end of sect. imi these Gaussian ensembles [3l|, may 
be realized by attaching an external harmonic spring, 
eq. ([H]), parallel to the system. As we have already seen 
in fig. [T1 this allows to gradually reduce the volume fluc¬ 
tuations while keeping constant the average pressure P 
of the system, i.e. the thermod yna mic (average) state of 
the system remains unchanged |34| . Expressed as a func¬ 
tion of A, the mean-squared volume fluctuations decrease 
linearly in agreement with eq. (1481) . Confirming eq. dsni), 
a similar scaling is also found for other moments of the 
volume (not shown). 

Focusing on simple ID nets the verification of several 
predicted pressure and volume-pressure correlation func¬ 
tions is presented in fig. |S1 The data presented in panel 
(a) has been obtained for systems with identical spring 
constants {6k = 0) at a low temperature T = 0.01. Vol¬ 
ume and pressure fluctuations are thus highly correlated 
in the NPT-ensemble, i.e. Creglp ~ Ij and the fluctua¬ 
tions of the excess pressure yp.exli in the NVT-ensemble 
is consequently small. The second system presented 
in panel (b) corresponds to a polydispersity 6k = 0.9 
and a temperature T = 1 for which Creglp ~ 0-44 and 
?7F,exli /??A ~ 0.8. Note that the vertical axes are made 
dimensionless by rescaling the data using either the com¬ 
pression modulus K or the affine dilatational elasticity 
ryA obtained for A = 0. We stress first that the com¬ 
pression modulus K may be computed irrespective of A 
by linear regression (spheres) of the measured {V, P) as 
predicted by eq. (l5^ . The corresponding dimensionless 
correlation coefficient Creg|;i, indicated by the squares if 
found to decrease monotonously from its maximum at 
A = 0 to zero in the NVT-ensemble. In agreement with 
eq. (PI) the decay is more sudden for values K/tja close 
to unity, i.e. if the non-affine contribution is negligible 
as in panel (a), and becomes more and more gradual for 
smaller K/tja as seen in panel (b). The linear decay of 
the total pressure correlation function rjp \ \ demonstrates 




FIG. 8: Verification of various predictions for generalized A- 
ensembles for simple ID nets: (a) identical springs {6k = 0) 
at temperature T = 0.01 and pressure P = 0 where K = 1, 
rjA = 1.02, Cregig « 1 and fo — 0.0199 for x — PidjK = 
0.01 and (b) polydisperse springs with Sk — 0.9 at T = 1.0 
and P = 0 where K = 0.6, rjA = 3.00, Creglg = 0.44 and 
fo = 0.53 for X = PidfK = 1.68. The lines compare the data 
with the respective prediction: the thin solid lines show the 
compression modulus obtained using the regression relation, 
eq. (1531) , the dash-dotted lines the transformation relation for 
the pressure fluctuations eq. (pl). the bold solid lines 

the key claim, eq. (©, for the Rowlinson functional Adrow|^. 


that the predicted transformation relation of the total 
pressure fluctuations, eq. dSl]), indeed holds as shown by 
the dash-dotted lines. Note that in the limit A —?> 1 
we have ?7 fIa ^F,exli due to the MC-gauge used for 
the instantaneous ideal pressure. Since for the first sys¬ 
tem in panel (a) the thermal noise is very week, we have 
w ?7F,ex|;^ for all A (not shown). Therefore, ryF.exI;^ 
is only presented in the second panel. Our key prediction 
eq. (0 for the Rowlinson functional is confirmed by the 
crosses presenting -pRow|;^/Ar vs. A. 

The scaling of ProwIa is further investigated in the 
main panel of hg. [5] for a broad range of ID nets (with 
parameters as indicated in the figure) and 2D pLJ beads 
(spheres) for one state point {T = 0.001, P = 2). Ac¬ 
cording to eq. 0) Jqrow|>, / AT — /o (a;) 0 in the NPT- 
limit. Reshuffling the contributions indicated in eq. 0, 
the latter limit has been used for the vertical axis. Using 
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FIG. 9: Scaling with A for various compression moduli K and 
reduced ideal pressures x — Pid/K. Main panel: Rescaled 
Rowlinson functional 7mow|;^ for different polydispersity Sk, 
pressure P and temperature T as indicated. The spheres in¬ 
dicate 2D pLJ beads, all other symbols data obtained for ID 
nets. Note that the horizontal scaling variable A(1 — fo{x)) = 
A(1 — x)'^ > 0. The bold line corresponds to eq. ([5|. Inset: 
Rescaled data for r?p,mix|;^ as a function of Xx{x — 1). The 
bold bisection-line corresponds to eq. (l58ll . Only data for ID 
nets with x > 1 are presented. 

as scaling variable A(1 — fo{x)) = A(1 — x)^ for the hor¬ 
izontal axis, all data must collapse on the bisection line 
according to eq. 0. As may be seen from the double- 
logarithmic plot, this is indeed the case for a broad range 
of data sets. The inset of fig. [S] presents a similar scaling 
plot for the function ?7F,mix|> characterizing the correla¬ 
tions of ideal pressure and excess pressure contributions 
for ID nets. The vertical axis corresponds again to the 
prediction of the NPT-limit, eq. ([35]). The horizontal 
axis Aa:(a: — 1) is chosen such that according to eq. (1551) 
all rescaled ?7F,mix|;^ have to fall on the bisection line. 
(Since the horizontal axis is logarithmic, only data with 
X > 1 can be represented and, hence, only data for ID 
nets are given.) This is confirmed by the presented data. 

V. CONCLUSION 

We have revisited in this paper theoretically and nu¬ 
merically various correlations of the normal pressure P 
and its contributions Pid and Pex for isotropic solids 
and fluids using simple coarse-grained models [l^, e.g., 
strictly ID networks of permanently fixed springs or the 
KA model for binary mixtures in three dimensions [3^. 
Making more precise several statements made in the ap¬ 
pendix of ref. (l3l | and extending the brief communication 
ref. [1^ , we have compared fluctuations in generalized A- 
ensembles where the volume fluctuations are tuned by 
means of an external harmonic spring potential allowing 
to switch gradually between the standard NPT- (A = 0) 
and NVT-ensembles (A = 1). We have stressed that the 
widely used stress fluctuation formula, eq. for the 


compression modulus K in the NVT-ensemble may be 
obtained directly 

• without the affine volume rescaling trick for the 
NVT-ensemble (reminded in appendix|B]) used first 
for liquids by Rowlinson Q and, 

• more importantly, without assuming a reference 
position for the particles and a microscopic dis¬ 
placement field which is only possible for solids 

using the general thermodynamic transformation rules 
between conjugated ensembles [2^ and assuming the sys¬ 
tems to be sufficiently large {V —>■ oo) for the given tem¬ 
perature T and compression modulus K . The direct 
thermodynamic derivation can readil y b e adapted to the 
shear modulus G in isotropic systems [Q and to the more 
general elastic moduli characterizing anisotropic solids as 
reminded in appendix |C] [2^ . 

The Rowlinson stress fluctuation functional -Frow [HI 
has been computed deliberately in the unusual NPT- 
ensemble (A = 0) to make manifest the general transform 
eq. (HI at the heart of the stress fluctuation formalism. 
We have demonstrated that FrowIq = ^’id(2—Pid/AT), i.e. 
FrowIo vanishes in the low-temperature limit. More gen¬ 
erally, we have investigated FrowIa other correlation 
functions as a function of the parameter A characterizing 
the volume fluctuations. As announced in the Introduc¬ 
tion, eq. the Rowlinson functional is found to interpo¬ 
late linearly between the classical ensemble limits (fig. [8]). 
Note that the specification of the compression modulus 
K and the dimensionless regression coefficient Creglg im¬ 
plies 77 a and 77F,exli (and vice versa) and together with 
the ideal pressure Pjd (which implies x = Pid/K) this 
allows the complete description of all the discussed cor¬ 
relation functions at different A. Our theoretical and nu¬ 
merically results, especially eq. ®, may allow to read¬ 
ily calibrate (correctness, convergence and precision) the 
various barostats commonly used pm. In the near 
future we plan (i) to also consider systems close to a 
first-order (e.g., solid to liquid) phase transition gener¬ 
alizing the work by Hetherington [Hl| and (ii) to extend 
our Mproach to negative values of A (following in that 
ref. [7|) increasing artificially the — then not necessar¬ 
ily Gaussian — fluctuations of the extensive variable and 
making thus the system increasingly unstable. 
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Appendix A: Affine interaction energy 


Affine displacement assumption. We consider here 
the change of the conservative interaction energy Us{V) 
of a configuration s under an imposed dilatational strain 
assuming that all particles respond to the macroscopic 
constraint by an affine microscopic displacement. It is 
convenient to introduce a dimensionless scalar character¬ 
izing the relative volume change 

e = 5V{e)/V{t)) = V{e)/Vffi)-l (Al) 


with V (0) being the instantaneous volume of the unper¬ 
turbed reference simulation box (e = 0). Obviously, for 
canonical ensembles the instantaneous volume V can be 
replaced by V. The affinity assumption implies, e.g., 

a;(0) ^x(e) =cc(0)(l + e)i/^ (A2) 


for the ^-coordinate of particle positions or relative par¬ 
ticle distances where the argument (0) denotes again the 
reference system. Equation (IA2I) implies that the squared 
distance between two particles transforms as 


r2(0) ^ r^[e) = r^(0)(l -f ef/^. 
It follows that 


dffi{e) 

de 

de^ 


2(2 - d) 


r^{0) 


(A3) 


(A4) 

(A5) 


where we have taken finally for both derivatives the limit 
e —)• 0. For the first two derivatives of a general function 
/(r(e)) with respect to e this implies 


df{r{e)) 

de 

d^fjrje)) 

de‘^ 




(A6) 

(A7) 


for e —)> 0 and dropping the argument (0) on the right 
hand-sides. 

Pair potential assumption. Let us assume now in ad¬ 
dition that the interaction energy is given by 


Us{e) =^u{ri{e)) (A8) 

i 


where the index I labels the interactions between the par¬ 
ticles i and j with i < j. We note en passant that for the 
first derivative of the interaction energy Us with respect 
to the volume eq. (IA6I) implies 

- U'siy) = -U's{e)/V{t)) = -^ E (^9) 


where we have taken e —>■ 0 in the last step. As stated in 
the main text, eq. ®, this is exactly the Kirkwood virial 
for the instantaneous excess pressure Pex- As defined 
in eq. ®, the coefficient r]A,ex/V measures the second 
derivative of the interaction energy Us{V) of a configura¬ 
tion s with respect to an affine dilatational strain. Using 
eq.dMI it follows for the instananeous value fjA.ex that 


VA.ex/V 


U’JiV) 


1 

U(0)2 


U's'ie) 


1 


E 


rfu"{ri) + riu{ri) 


1 


'^nu'iri) 

I 


(AlO) 


where we have taken again e —>■ 0 and have dropped the 
index (0) in the last step. 

Simple averages. Only the affinity assumption and 
the pair potential choice have been used up to now. The 
mean value r/A.ex = (^A,ex) is then obtained by taking 
the thermal average over all configurations s and over all 
volumes V depending on the ensemble. Since all con¬ 
tributions to eq. (lAIOl) correspond to simple averages, 
thermostatistics tells us that one can replace V for suffi¬ 
ciently large systems by its mean value V, eq. m- One 
confirms that the affine dilatational elasticity becomes 


^A,ex 


^Born 


??Born + Pex with (All) 

(^.rfu"{ri) + riu\ri) \ (A12) 


as already stated in eq. O- Note that it is inconsistent to 
neglect the explicit excess pressure contribution to r]A,ex 
in eq. (IA11|) but to keep the underlined contribution to 
??Born which amounts to —Pex./d. The sum of both terms 
Pexffi — l)/d only vanishes in d = 1 dimensions where 


VA,ex — Y 



as used in sect. IIVBI 


(A13) 


Appendix B: Volume rescaling trick revisited 

Introduction. The stress fluctuation formula, eq. ®, 
for the compression modulus K of isotropic systems has 
been derived in sec. Ill Cl essentially using the general 
transformation relation eq. between conjugated en¬ 
sembles and evaluating the pressure fluctuations in the 
NPT-ensemble by integration by parts, eq. (1^^ . Histori¬ 
cally, eq. m has been first derived properly by Rowlinson 
Q who computed using the volume rescaling trick the 
second derivative, eq. o, of the free energy F{T, V) = 
—kBT\og[Z[V)\. (The NVT-ensemble is thus used and, 
hence, U = U in the following.) Since the total system 
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Hamiltonian may be written as the sum of a kinetic en¬ 
ergy and a potential energy, the partition function Z{V) 
factorizes in an ideal contribution Zi^{V) and an excess 
contribution Zex(H) on which we focus below. We re¬ 
mind 0 that using Zi^{V) ^ one readily confirms 
^id = k^Tp and = ^id for the ideal contributions 
to, respectively, the total pressure P = Pid + Pex and the 
total compression modulus K = Ki^ + K^x- 

Mapping of strained and unstrained eonfigurations. 
Since both the perturbed as the unpertubed partition 
function is sum over all possible particle configurations, 
one can always compare the contribution of a configu¬ 
ration s of the strained system with the contribution of 
a configuration of the reference being obtained by the 
affine rescaling of all coordinates according to eq. (IA2I) . 
Using the same notations as in appendix the interac¬ 
tion energy Us{e) of the strained system can be expressed 
in terms of the coordinates (state) of the unperturbed 
system and the explicit metric parameter e. The ex¬ 
cess contribution Pex(e) = —fcBT'-^ex(e) is thus obtained 
from Zex(e) = Ss where constant prefac¬ 

tors have been omitted and where the sum is taken over 
all possible configurations s. 

General conservative potential. We note for the first 
two derivatives of the free energy 


[ln(Zex(e))]' = 
[ln(Zex(e))]" = 


KM 

■^ex(e) 


2 


(Bl) 

(B2) 


and of the partition function 


Being a simple average the first term on the right hand- 
side corresponds to r?A.ex as defined in eq. dSD). Using 
eq. (1551) the second term is seen to yield the pressure 
fluctuation contribution ?7F,exli- Summarizing all terms 
we have thus shown that 


K = Kid + Kex = Pid + »?A.ex — ''fF.exli (B7) 

in agreement with eq. m if the Born-Lame coefficient is 
defined more generally as rjBom = ??A.ex — Pex- 

Pair potential choice. Up to now we have stated Pex 
and PTex for a general interaction potentiel Us{e). Assum¬ 
ing the interactions to be described by a pairwise additive 
potential, eq. (IA8I) . the Kirkwood relation, eq. , is con¬ 
firmed using eq. (IB5I) and eq. (IA9|) . As already noted at 
the end of appendixone confirms using eq. (lAlOl) that 
??Born = T?A.ex “ Pex agrees with eq. (g]). 


Appendix C: General stress fluctuation formalism 

Introduction. We remind here the general stress fluc¬ 
tuation formalism derived by Squire, Hold and Hoover Q 
and show that Rowlinson’s formula, eq. ®, is a special 
case obtained by symmetry considerations. For conve¬ 
nience we introduce the two linear projection operators 

1 ^ 

r^[A^P] = (Cl) 

1 

1 ^ 

^ Y (C2) 

a,^, 7 ,(5=1 


zM) 

z'M 


S 

YiPUMf 

s 

S 


(B3) 


(B4) 


with a prime denoting again a derivative with respect to 
the indicated argument. Using Pex = —dFex{V)/dV and 
eq. (IB3I) and taking the limit e —>■ 0 one verifies that 

Pex = (-Pex^ with Pex = UM\e=0 (B5) 

which defines the instantaneous excess pressure. (The av¬ 
erage taken uses the weights of the unperturbed system.) 
The excess pressure thus measures the average change of 
the total interaction energy Us{e) taken at e = 0. The 
excess compression modulus Kg^ is obtained using in ad¬ 
dition eq. (IB2I) and eq. (IB4I) and taking finally the e —>■ 0 
limit. This yields 


pVKg^ 




(B6) 


with A°‘^ and being, respectively, second- and 

forth-rang tensors and dap the Kronecker symbol [s^. 
Greek letters are used for the spatial coordinates 
a,l3,j,d = l,...,(i. The following identities are read¬ 
ily verified 

= pM] r2[p“^] (C3) 

r4A°‘M = 1 for (C4) 

r4A°‘M = l/diov = Sa,.ySps (C5) 
r4A°‘M = l/d^ for = (C6) 

with n“,... being the spatial components of a normal¬ 
ized vector, i.e. = 1- 

Thermodynamics and symmetry. Generalizing the 
definition of the pressure P given in eq. o, the stress 
tensor cr“^ may be defined as the first derivative of the 
free energy per volume with respect to the linear strain 
e“^ [l[ characterizing the macroscopic deformation of the 
system. Note that in general does not vanish at the 
investigated state point, i.e. the systems may be pre¬ 
stressed at the reference strain = 0. The latter point 
is crucial for essentially all soft matter systems which 
only assemble because a finite density and/or stress is 
applied. It is less important for the classical crystalline 








15 


solids where the elastic moduli are normally huge com¬ 
pared to the imposed stresses. We denote by 5a°‘^ the 
increment of the stress tensor cr“^ under an increment 
generalizing the dilatational strain e used in ap¬ 
pendix [BJ Assuming an infinitessimal strain increment, 
Hooke’s law reads quite generally [l| 

d 

^£7-5 (C7) 

7,(5—1 

where the elastic moduli i<;“^7<5 stand for the second 
derivative d^F/de°‘^de^^ of the free energy per volume at 
the given thermodynamic state [T| . Let us now assume a 
pure dilatational strain without shear, i.e. = e Sap. 
As may be seen from eq. (4.6) of ref. [Ij, this implies 

= dKeSap. Hence, (5cr““ = cPKe. Or, using 

eq. (IC7I) one sees that = J2ap Com¬ 
paring both expressions, this shows that the compression 

modulus is given by AT = 74 using the linear pro¬ 

jection operator defined above. 

Stress fluctuation relations. As described in the liter¬ 
ature i i m, na H, the elasticity tensor can 

numerically be computed from the sum 

^ +^Born+^F.id + ^F.ex 

with contributions as specified below. The compression 
modulus is obtained by applying the linear operator Ti 
to each term and by summing up the contributions. Note 
that some authors only refer to the sum + C‘^avn + 

^F^d + ^F ex^ elasticity tensor [1, d, [T^l ■ 

Initial stress contribution. The (often not included) 
first term _B“^7<5 gq^ (|C8|) may we written as [T^ 

^a/375 = _^o.l3s^S+Cr°‘'^Sps+<y°‘^Sp.y (C9) 

= P [SapSjS — SajSps — SasSpj) . (CIO) 

Following Birch we have assumed in the second step 
that the system is isotropically stressed, i.e. cr“^ = 
—PSap. Since 5“^^^ becomes negligible for small to¬ 
tal pressures P ~ 0, this term is often not computed. 
Consistency implies then, however, to set Pid = —Pex in 
the remaining terms of eq. (IC8I) . Returning to eq. (ICIOI) 
one sees using eq. (El and eq. (IC5I) that 

7^j^a/375l (Cll) 

\ dj d d 

There is thus no contribution to K from the Birch term 
for d = 2 dimensions and a contribution — P for d = 
1. The leading term P in the second step of eq. (ICllI) 
corresponds to the total pressure indicated in eq. (0. 

Ideal pressure contributions. The second and the 
forth term in eq. (IC8I) correspond to the kinetic (ideal) 
contributions 

^a/375 = 2Pii{Saj6ps + SasSp^), (C12) 

= -P,^{5a^Sps + Sa5Sp^). (C13) 


Using again eq. (IC4I) and eq. (IC5I) this yields T 4 = 

4Pid/d and Ti= — 2Pid/d. Together with the 
ideal pressure contribution from eq. (ICllI) all the kinetic 
terms sum up (independently of the dimension) to Pid in 
agreement with eq. ©• It is worthwhile to stress that if 
the Birch term P“^75 jg neglected for Pjd —Pexj the 
ideal pressure contributions sum up incorrectly as is the 
case in ref. Q. In fact, the traditional separation of the 
(rather trivial) different kinetic terms stated in eq. (IClOll 
is from the numerical point of view unfortunate. More 
importantly, this separation is theoretically misleading: 
The total system Hamiltonian being the sum of the ki¬ 
netic energy and the conservative interaction potential, 
all kinetic contributions to the elastic moduli can (and 
should) be factorized out from the start. This leads im¬ 
mediately to the obvious total kinetic contribution 

= P.ASapS^p (C14) 

replacing the aforementioned three terms. 

Affine Born contribution. Generalizing the Born co¬ 
efficient 77Born in eq. © , the third term in eq. (IC8II 

corresponds to the (second derivative) of the energy 
change assuming the affine displacement of all particles 
(as in appendix lAl for a pure dilatational strain). Assum¬ 
ing pair interactions it is given by 

CsorTf = ^ ~ P«'(P)) 

(CI5) 

where n“, ... stand for the normalized vector between the 
pair of interacting particles labeled by 1. Using eq. (IC6I) 
one confirms that 

= ■^(^rfu'\ri)-riu'{ri)^ 

= 7Born + (C16) 

a 

where we have used the Kirkwood formula, eq. ®, in the 
last step. Note that the underlined term cancels exactly 
the underlined last term in eq. (ICllI) . This is equivalent 
to the statement made after eq. (IA12I) in appendix lAl 
Stress fluctuations. The last two terms indicated in 
eq. (ICSp stand for the correlation function 

C^PjS ^ ^^PyS ^ ^a/375 = _py (C17) 

with denoting the fluctuation of the total stress ten¬ 
sor = a^f +cr“/. Using the factorization of ideal and 
excess stress fluctuations in the canonical ensemble at 
constant volume V (as shown in sect. Ill Ell this is different 
in NPT-ensembles), both contributions can be separated. 
The ideal stress fluctuation contribution has already been 
accounted for above. We remind that the instantaneous 
excess pressure Pex is given by the trace Pex = T^iPex^] 
over the excess pressure tensor Pg“^ being the negative 
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of the instantaneous excess stress tensor. Using eq. (IC3I) 
one thus obtains 

= -pv {r^[5PsiY) (C18) 

which is identical to the pressure fluctuation term in 
eq. ©. We have thus confirmed that the general stress 


fluctuation relation, eq. (IC8I) . is consistent with the Row! 
inson formula for the compression modulus K. 
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